Section 1- Preliminaries

Load libraries

To begin with the analysis, we will first load the required packages. The downstream analysis of the scRNA-seq is built upon Seurat and sctransform packages. The scDblFinder package is used for identifying and removing doublets in the dataset. The single-cell annotation is carried out using the marker-based scSorter method.

library("Seurat")
library("ggplot2")
library("sctransform")
library("BiocParallel")
library("scDblFinder")
library("SingleCellExperiment")
library("cowplot")
theme_set(theme_cowplot())
library("RColorBrewer")
library("viridis")
library("scSorter")
library("knitr")

Load data

First, we will load the filtered count matrix from cellranger. This matrix contains the scRNA-seq counts from all the HPAP samples.

panc_data <- Read10X(data.dir = "./data/HPAP_samples_aggr_18Apr2022/outs/count/filtered_feature_bc_matrix/")

Dimension of raw merged counts matrix

The dimension of the raw counts matrix after merging all HPAP scRNA-seq samples with cellranger is shown below. There are total of 36601 genes and 269417 cells.

# panc_data <- Read10X(data.dir = "PATH_TO_FEATURE_MATRIX")
dim(panc_data)
## [1]  36601 269417

Section 2- Preprocessing the data

Create Seurat object

The Seurat object initialization step above only considered cells that expressed at least 200 genes. The min.cells- parameter will include features detected in at least three cells. The min.features - Include cells where at least three fifty features are detected.

With the above parameters the count matrix will also be subsetted. Optionally, to revive the excluded features, create a new object with a lower threshold.

panc <- CreateSeuratObject(counts = panc_data, min.cells = 3, min.features = 200,
                           project = "hpap@67")
dim(panc)
## [1]  32066 258379

In the Seurat object, 32066 genes and 258379 cells were present after passing the thresholds.

Add metadata

Here, we will add four new columns to the Seurat metadata. The first column sample_id contains raw sample names from cellranger, hpap_id refers to the donor identifier, the disease_id corresponds to the numeric count in each disease type, and finally the disease_state shows the general state of disease (CTL, AAB, and T1D).

## Add metadata information
sample <- read.csv(file.path("./metadata", "HPAP_samples_aggr_18Apr2022_cellranger.csv"), stringsAsFactors=F)
## Add sample_id
CellsMeta["sample_id"] <- df$sample_id
CellsMetaTrim <- subset(CellsMeta, select = c("sample_id"))
panc <- AddMetaData(panc, CellsMetaTrim)

## Add hpap_id
CellsMeta["hpap_id"] <- df$hpap_id
CellsMetaTrim <- subset(CellsMeta, select = c("hpap_id"))
panc <- AddMetaData(panc, CellsMetaTrim)

## Add disease_id
CellsMeta["disease_id"] <- df$disease_id
CellsMetaTrim <- subset(CellsMeta, select = c("disease_id"))
panc <- AddMetaData(panc, CellsMetaTrim)

## Add disease_state
CellsMeta["disease_state"] <- df$disease_state
CellsMetaTrim <- subset(CellsMeta, select = c("disease_state"))
panc <- AddMetaData(panc, CellsMetaTrim)

Table showing metadata information

head(panc@meta.data) %>%
  kable(format = "html", col.names = colnames(head(panc@meta.data))) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "100px")
orig.ident nCount_RNA nFeature_RNA sample_id hpap_id disease_id disease_state
AAACCTGAGCGTTGCC-1 16875 4080 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D
AAACCTGAGGACAGAA-1 6497 892 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D
AAACCTGAGTCGTACT-1 1169 253 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D
AAACCTGCAACTGCTA-1 13689 3553 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D
AAACCTGCAATGTAAG-1 16375 3583 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D
AAACCTGCACCAGGTC-1 16431 3438 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D
## Calculate MT %
panc[["percent.mt"]] <- PercentageFeatureSet(panc, pattern = "^MT-")

The below plot shows the number of features, counts, and mitochondrial reads before filtering.

VlnPlot(panc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
        pt.size = 0.001, combine = T) & theme(legend.position = 'none',
                                                               axis.title.x = element_blank())

VlnPlot(panc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
        pt.size = 0, combine = T) & theme(legend.position = 'none', 
                                                           axis.title.x = element_blank())

Correlation between nCount_RNA - nFeature_RNA and nCount_RNA - percent.mt before filtering.

# FeatureScatter used to visualize feature-feature relationships.
plot1 <- FeatureScatter(panc, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "disease_state")  & theme(legend.position = 'none')
plot2 <- FeatureScatter(panc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "disease_state") 
CombinePlots(plots = list(plot1, plot2))

Remove doublets

The scDblFinder function is employed to remove the potential doublets from each HPAP donor separately. The parameter sample_id refers to individual donor. The function uses an object from SingleCellExperiment class, therefore the SingleCellExperiment package is used for conversion.

## Convert object into singlecellexperiment
panc.sce <- as.SingleCellExperiment(panc)
panc.sce <- scDblFinder(panc.sce, samples="sample_id", clusters=FALSE, BPPARAM=MulticoreParam(10))
## Convert sce object back to seurat
panc_seurat <- as.Seurat(panc.sce, counts = "counts", data = "logcounts")

The total number of doublets are shown below

table(panc.sce$scDblFinder.class)
## 
## singlet doublet 
##  234285   24094

The doublet information is stored in the metadata. The updated metadata can be found below

Updated metadata table with doublets information

head(panc@meta.data) %>%
  kable(format = "html", col.names = colnames(head(panc@meta.data))) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "100px")
orig.ident nCount_RNA nFeature_RNA sample_id hpap_id disease_id disease_state percent.mt
AAACCTGAGCGTTGCC-1 16875 4080 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D 7.330370
AAACCTGAGGACAGAA-1 6497 892 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D 1.062029
AAACCTGAGTCGTACT-1 1169 253 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D 4.362703
AAACCTGCAACTGCTA-1 13689 3553 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D 3.382278
AAACCTGCAATGTAAG-1 16375 3583 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D 3.877863
AAACCTGCACCAGGTC-1 16431 3438 HPAP021_75162_T1DIslets HPAP021 T1D_ 1 T1D 1.545858

Pre-processing the data based on information from above steps

The nFeature_RNA and nCount_RNA represents the number of genes detected in each cell and the total number of molecules (UMIs) detected within a cell respectively. While the low nFeature_RNA in a cell mean that it may be dead/dying or it may represent an empty droplet, the high nCount_RNA and nFeature_RNA denotes that the cell may be a doublet or multiplet. This filtering along with mitochondrial reads is crucial pre-processing step, because, removing such outliers from these groups might also remove some of the doublets or dead/empty droplets.

Remove all the potential doublets found

The cell being either doublet or singlet is stored in column scDlFinder.class. The potential doublets are removed in this step.

# After doublet removal
panc <- subset(panc_seurat, subset = scDblFinder.class  == "singlet")

The data after filtering based on nFeature and percent.mt is shown below

panc_new <- subset(panc, subset = nFeature_RNA > 200 & nFeature_RNA < 9000
                   & percent.mt < 25
                   & nCount_RNA < 120000
)
dim(panc_new)
## [1]  32066 222679

Update Metadata with MT and Rb genes information

After removing doublets, now we can proceed with further pre-processing steps. For ex, Mitochondrial genes are useful indicators of cell state.

Add the mitochondrial proteins information to the metadata after fitering

Adding the Mitochondrial genes

panc[["percent.mt"]] <- PercentageFeatureSet(panc, pattern = "^MT-")

Display the metadata information after filtering

The below plot shows the number of features, counts, and the percentage of mitochondrial reads after filtering.

VlnPlot(panc_new, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
        pt.size = 0.001, combine = T) & theme(legend.position = 'none',
                                                               axis.title.x = element_blank())

VlnPlot(panc_new, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
        pt.size = 0, combine = T) & theme(legend.position = 'none', 
                                                           axis.title.x = element_blank())

Correlation between nCount_RNA - nFeature_RNA and nCount_RNA - percent.mt after filtering.

The FeatureScatter is typically used to visualize feature-feature relationships.

plot1 <- FeatureScatter(panc_new, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "disease_state")  & theme(legend.position = 'none')
plot2 <- FeatureScatter(panc_new, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "disease_state") 
CombinePlots(plots = list(plot1, plot2))

Section 3- Normalizing and scaling the data

SCT: Single Cell Transform

The seurat’s SCTranform (SCT) function is used to normalize the counts that measures the differences in sequencing depth per cell for each sample. The SCT method is built on the concept regularized negative binomial model to perform normalization and variance stabilization of single-cell RNA-seq data. It removes the variation due to sequencing depth (nUMIs). The vars.to.regress parameter allows us to regress out the variation from other sources such as percentage of mitochondrial. The output of SCT model is the normalized expression levels for all the transcripts. Lastly, the variable.features.n parameter is used to select the variable features and is set as 3000.

panc_new_sct <- SCTransform(panc_new, method = "glmGamPoi", variable.features.n = 3000, vars.to.regress = "percent.mt", verbose = FALSE)

The seurat object now consists of two assays RNA and SCT, with the default being set as SCT.

DefaultAssay(panc_new_sct) 
## [1] "SCT"

Perform dimensionality reduction by PCA and UMAP embedding

These are now standard steps in the Seurat workflow for visualization and clustering

panc_new_sct <- RunPCA(panc_new_sct, verbose = TRUE)

Plot showing heatmaps of PCs- 1 to 12 from the PCA analysis

DimHeatmap(panc_new_sct, dims = 1:12, cells = 500, balanced = TRUE)

Plot showing heatmaps of PCs- 13 to 24 from the PCA analysis

DimHeatmap(panc_new_sct, dims = 13:24, cells = 500, balanced = TRUE)

Plot showing heatmaps of PCs- 25 to 36 from the PCA analysis

DimHeatmap(panc_new_sct, dims = 25:36, cells = 500, balanced = TRUE)

Plot showing the variance explained by PCs from the PCA analysis

ElbowPlot(panc_new_sct, ndims = 40)
ElbowPlot(panc_new_sct)

Run UMAP

The Uniform Manifold Approximation and Projection (UMAP) dimensional reduction technique is run on the first 20 principal components.

panc_new_sct <- RunUMAP(panc_new_sct, dims = 1:20, verbose = FALSE)

Find Neighbors

Here we compute the k nearest neighbors in the dataset. The nearest neighbors are calculated on the reduced dimensions from the PCA. We use first 20 principal components.

panc_new_sct <- FindNeighbors(panc_new_sct, dims = 1:20, verbose = FALSE)

Find Clusters

Once we have calculated the k-nearest neighbors in previous step, we will find the clusters of cells by the shared nearest neihbor (SNN) clustering algorithm. Here, we have set the resolution as 0.8. The higher resolution usually gives more number of clusters.

panc_new_sct <- FindClusters(panc_new_sct, resolution = 0.8, method = "igraph", verbose = FALSE)

UMAP plotting

UMAP plot showing the clusters based on SCT data and snn optimization at resolution 0.8

DimPlot(panc_new_sct, group.by='SCT_snn_res.0.8', reduction='umap') +
  ggtitle('SCT_snn_res.0.8')

UMAP plot showing the distribution of cells across different disease types.

DimPlot(panc_new_sct, group.by='disease_state', reduction='umap') +
  ggtitle('Disease type')

UMAP plot showing the cell distribution with respect to individual donors

DimPlot(panc_new_sct, group.by='hpap_id', reduction='umap') +
  ggtitle('Sample type') 

Section 4- Cell Annotation

scSorter: single cell Sorter

Loading marker genes data

anno <- read.csv("./annot/anno.csv")
anno[,1] <- NULL

Display marker genes across celltypes

We used marker-based scSorter cell annotation method for assigning cells to known cell type based on the marker genes. The markers for each cell type were added based on comprehensive literature search.

anno %>%
  kable(format = "html", col.names = colnames(anno)) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "100px")
Type Marker
Acinar PRSS1
Acinar REG1A
Acinar CPA1
Acinar CPA2
Alpha GCG
Alpha GC
Alpha TTR
Beta INS
Beta IGF2
Beta IAPP
Ductal KRT19
Ductal CFTR
Ductal SFRP5
Ductal MMP7
Stellates_Mesenchymal COL1A1
Stellates_Mesenchymal PDGFRB
Stellates_Mesenchymal RGS10
Stellates_Mesenchymal THY1
Endothelial VWF
Endothelial CD93
Immune NCF2
Immune PTPRC
PP_Gamma PPY
PP_Gamma CARTPT
PP_Gamma PCDH10
PP_Gamma PLAC8
Delta SST
Delta LEPR
Delta PRG4
Epsilon GHRL

The scSorter method is based on a semi-supervised learning algorithm. It expects the pre-processed data as input. Therefore, we used the previously normalized, scaled and transformed expression matrix generated by SCT method as input. The top 3000 highly variable genes selected by SCT method was also given as input for scSorter method.

expr <- panc_new_sct@assays$SCT@data
topgenes <- head(panc_new_sct@assays$SCT@var.features, 3000)
topgene_filter = rowSums(expr[topgenes, ]!=0) > ncol(expr)*.1
topgenes = topgenes[topgene_filter]

## At last, we subset the preprocessed expression data and run scSorter.
picked_genes = unique(c(anno$Marker, topgenes))
expr = expr[rownames(expr) %in% picked_genes, ]

fit the scSorter method

Based on the expression matrix and the marker genes the scSorter method is fit and the cell type assignment results predicted are stored in the Pred_Type vector.

rts <- scSorter(expr, anno)
panc_new_sct$cell_type <- rts$Pred_Type
table(rts$Pred_Type)
## 
##                Acinar                 Alpha                  Beta 
##                 59496                 65075                 45009 
##                 Delta                Ductal           Endothelial 
##                  5664                 20300                  5864 
##               Epsilon                Immune              PP_Gamma 
##                   196                  1220                  2821 
## Stellates_Mesenchymal               Unknown 
##                  9583                  7451

Add celltype information into the metadata

The celltype information is added to the metadata in the original Seurat object

cell_type.info <- data.frame(cell_type = panc_new_sct$cell_type, row.names= colnames(panc_new_sct_umap))
panc_new_sct_umap <- AddMetaData(object = panc_new_sct_umap, metadata = cell_type.info)

## Convert cell type to factor
panc_new_sct_umap$cell_type <- as.factor(panc_new_sct_umap$cell_type) 

Idents(panc_new_sct_umap) <- "cell_type"
panc_new_sct_umap$ident <- NULL

Check expression of markers

The expression of markers across all the cell types are shown in different UMAPs

FeaturePlot(panc_new_sct_umap, features = c("PRSS1", "REG1A", "CPA1", "CPA2", "GCG", 
                                             "GC", "TTR", "INS", "IAPP", "GHRL",
                                             "COL1A1", "PDGFRB", "RGS10", "THY1", "VWF",
                                             "CD93", "PPY", "SST", "NCF2", "PTPRC"), pt.size = 0.2, ncol = 5) + theme(legend.position = 'none', 
                        axis.title.x = element_blank(), axis.title.y = element_blank())

The UMAP showing the cell type annotations from the scSorter method

DimPlot(panc_new_sct_umap, group.by = "cell_type", reduction = "umap", label = FALSE, cols= c("#e30800", "#f56505", "#dec400", "#006630", "#0223c7","#5b02c7", "#00b0e6", "#c40080", "#02f00a", "#7d3301", "#000000"))

Below we show the UMAPs for the top 5 cell types

UMAP showing only Acinar cells

DimPlot(panc_new_sct_umap, group.by = "cell_type", reduction = "umap", label = FALSE, cols= c("#fac720", "#000000", "#000000", "#000000", "#000000","#000000", "#000000", "#000000", "#000000", "#000000", "#000000", "#000000"))

UMAP showing only Alpha cells

UMAP showing only Beta cells

UMAP showing only Ductal cells

UMAP showing only Immune cells